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20  Abstract 

to-thickness  ratios  100  and  200  , respectively.  The  shapes 
of  the  buckled  shells  are  computed  for  the  first  time  in  an 
unexpected  hint  that  there  is  an  asymptotic  (ratios-to- thickness 
ration- independent)  buckled  shape.  Numerical  solution  is  by 
means  of  spectral  (Galerkin)  expansions  of  up  to  60  modes  in 
associated  Legendre  functions  of  order  one  applied  to  the 
quadratically  non-linear  version  of  E.  Reissner's  equations. 


CAPTIONS  TOU  FICUUKS 


Fig.  2 Graph  of  dimensionless  pressure  p versus  inward  polar 

deflection  (in  radii)  for  a/h  = 100,  both  symmetric  and  asymmetric 

branches  (see  text).  Computed  with  40  modes  through  the  deflection 
* 

corresponding  to  P , then  the  remainder  with  24  modes. 

Li 

Fig.  3 Analogue  of  Fig.  2 for  a/h  = 200.  Analogous  computations 
with  60  modes  and  40  modes  respectively . 

Fig.  4 Shape  of  planar  cross-section  (through  axis  of  symmetry) 

, * 

of  asymmetrically  deformed  shell  for  both  a/h  = 100  and  200  and  p 
at  the  respective  asymmetric  minima.  The  dashed  piece  of  circle  shows 

the  remainder  of  the  undeformed  shell;  the  shapes  of  the  symmetrically 

* 

deformed  shells  for  p = p are  not  shown  because  each  intersects 

L 

itself  and,  hence,  is  not  physically  realizable  a priori. 
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BUCKLING  OF  A COMPLETE  SPHERICAL  SHELL  UNDER  UNIFORM  EXTERNAL  PRESSURE 


"I 

! 

[i 

by 

Harry  E.  Itaucli,  Neal  H.  Jacobs,  and  Jonathan  L.  Mans 
1 . Introduction 

A brief  history  and  discussion  of  the  title  topic  appears  in  the 

I 

preliminary  sketch  of  this  research,  [ll],  where  reference  is  made  to 
ri] , [3],  [5],  [6],  [13],  and  [14],  The  present  paper  is  self-contained. 

The  immediate  goal  of  the  research  and  the  paper  is  to  compute  the 
lower  critical  pressure  of  a complete  spherical  shell  for  plausible  values 
of  radius-thickness  ratios.  We  deal  with  axi symmetric  buckling  only  here, 
and  this  assumption  is  retained  throughout  the  paper  without  further 
mention.  It  is  recalled  that  the  lower  critical  pressure  is  the  smallest 
(greatest  lower  bound) of  the  pressures  for  which  the  initially  perfect 
shell  assumes  an  axi symmetric  buckled,  i.e.,  non-spherical  shape  and  thus 
represents  a theoretical  absolute  least  failure  load  for  the  shell. under 
axisymmetric  deformation.  In  view  of  the  marked  imperfection  sensitivity 
of  the  shell  ([5],  [6],  [13],  and  Section  7 below),  the  lower  critical 
pressure  may,  under  certain  circumstances,  be  the  only  reliable  theoretical 
failure  load,  as  was  suggested  by  von  Karman  and  Tsien  in  their  pioneer 
work.  At  any  rate  the  low  values  obtained  here,  roughly  7%  and  5%  of 
the  classical  linear  buckling  loads  for  the  radius-thickness  ratios  100 
and  200,  respectively,  make  the  intuitively  sensed  strength  of  the 
complete  spherical  shell  seem  illusory.  Coupled  with  related  results 
for  the  axial  compression  of  circular  cylinders  ([7],  [4]),  this  in- 
dicates the  caution  necessary  in  the  uso  of  these  highly  symmetrical 
thin-walled  structures  in  contrast  to  flat  plates  under  edge  compression, 

which  exhibit  postbuckling  stability.  Indeed,  tho  discovery  that  certain 
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optimized  structures,  o.g.,  flat  panols  reinforced  with  stringors, 
exhibit  imperfoetion-sonsi tivity  when  undor  end  compression  indicates 
that  it  may  be  necessary  more  often  to  perform  non-linear  lower  critical 
load  analyses  in  addition  to  the  now  customary  upper  critical  load, 
i.e.,  bifurcation  or  snap-through  analyses. 

The  basic  equations  used  in  the  present  analysis  are  the  small 
finite  deflection  form  of  E.  Reissner's  coupled  pair  of  non-linear 


ordinary  differential  equations  for  axisymmetric  deformation  of  axi- 
symmetric  shells  [ 12] . 


The  method  of  solution  is  a version  of  the  Galerkin  or  spectral 
method  in  which  the  complete  set  of  functions  used  is  that  of  the  as- 
sociated Legondre  functions  of  order  one.  These  functions  are  the  formal 
eigenfunctions  of  the  differential  operator  appearing  in  the  differential 
equations  so  that  the  linearized  system  is  diagonal.  In  the  opinion  of 
the  authors  the  spectral  method  has  three  features  to  recommend  it: 

(a)  it  is  conceptually  simple  and  relatively  straightforward  to  apply, 

(b)  the  resulting  Fourier-type  analysis  of  the  relevant  functions  into 
the  various  modes  is  enlightening  and  relates  the  non-linear  analysis 
directly  to  a familiar  method  of  linear  analysis,  and  (c)  the  method 
extends  directly  to  certain  more  complex  situations,  where  partial 
rather  than  ordinary  differential  equations  govern. 

The  progress  to  be  reported  here,  above  and  beyond  that  in  Til], 
is,  first,  the  exhibition  of  formulas  in  closed  form  for  the  cubic 
integrals  (-t,m,n),  see  Section  3,  and  the  consequent  evaluation  of  them 
by  the  program  TABLE  (see  Section  4),  and,  second,  the  numerical  solution 
of  tho  coupled  quadratic  equations  obtained  from  tho  spectral  method  by 


means  of  the  program,  SPHERE  (sco  Section  4) . 


I . 


-3- 

As  mentioned  in  [ll],  tho  first-named  writer  had  pursued  the  present 
research  as  far  as  represented  there  independently,  at  which  time  he 
became  aware  of  [ l] , in  which  tho  same  differential  equations  (with 
different  interpretations  for  the  dependent  variables)  are  solved  by 
quite  a different  method,  parallel  shooting,  which  is  confined  in  principle 
to  ordinary  differential  equations.  The  choice  of  data  presented  here 
is  motivated  by  tho  desire  to  give  results  both  of  autonomous  interest 
and  of  sufficiently  different  character  from  those  in  [ l]  to  justify  the 
presentation  of  tho  present  method  as  an  alternative  (see  Section  5 for 
comparison  of  results).  Particular  attention  is  called  to  Figure  4 and 
the  relevant  discussion  in  Section  5 and  the  two  other  "experimental 
discoveries"  there. 

Some  remarks  on  a method  of  incorporating  imperfection-sensitivity 
into  the  theory  given  here  are  given  in  Section  7. 

2.  The  Basic  Equations 

The  source  is  Reissner  [12,  Sections  2-4,  9-11,  in  particular 
Eqs.  (63)-(69)].  Tho  middle  surface  of  the  undeformed  spherical  shell 
of  radius  a and  thickness  h is  represented  in  cylindrical  coordinates 
(r,z,9)  by  r = rQ  = a sin  §,  z = zQ  = -a  cos  where  | is  the  co- 

latitude as  shown  in  Fig.  1,  where  cross  sections  (say,  9=0)  of  the 
undeformed  and  deformed  shells  are  shown.  Assuming  axisymmetric  de- 
formation with  the  z axis  as  axis  of  symmetry,  one  can  represent  the 
middle  surface  of  the  deformed  shell  by  r = rQ  + u,  z = zQ  + w,  whore 
u and  w,  functions  of  §,  are  respectively  the  radial  (horizontal)  and 
axial  (vortical)  components  of  the  displacement  vector  (u,w).  It  is 
important  to  note  that  u < 0 represents  displacement  inward  toward 


i 
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the  axis  of  symmetry,  while  w >0  in  the  southern  hemisphere  and  w < 0 
in  the  northern  hemisphere  represents  displacement  inward  toward  the 
equatorial  plane  z = 0 . The  angle  between  the  radial  (horizontal) 
direction  and  the  ray  tangent  to  the  meridian  of  the  undeformed  middle 
surface  in  the  direction  of  increasing  § at  any  point  with  colatitude 
§ is  denoted  by  and  the  corresponding  angle  at  the  displaced  points 

(coming  from  those  with  eolatitude  §)  on  the  deformed  middle  surface  is 
denoted  by  tp.  (See  Fig.  1.) 

The  first  basic  dependent  variable  3 is  defined  by 

(1)  P = -up  - «?0>  • 

It  is  important  to  observe  that  in  the  situation  shown  in  Fig.  1,  i.e.,  an 
inward  directed  dimple  at.  the  south  pole  (north  pole)  , one  has  3 >0  (<  0)  . 
One  defines  the  stress  function  \ji  by 

(2)  ijr  = r0«  = a (sin  §)H  , 

where  II  is  the  horizontal  (radial)  stress  resultant  at  all  points  on  the 
deformed  shell  which  were  originally  specified  by  5 on  the  undeformed 
shell.  It  is  important  to  observe  that  near  the  south  pole  ijr  < 0 implies 
compressive  stress,  with  the  same  implication  for  the  opposite  inequality 
at  the  north  pole. 

If  the  shell  is  subjected  to  a uniform  inward  normal  pressure  p and 
is  in  the  membrane  state,  3 = 0,  then  one  has 

'll  = Pa^  cos  5 sin  ? . 

Now  one  defines  the  second  basic  dependent  variable  iji  by 

(3)  <1/  = \ pa  cos  § sin  § + 1|i 

so  that  iji  doscribos  the  deviation  of  the  stress  function  i|i  from  the 
membrane  state. 


14 


Now,  if  D = Eir/12(l  - v ) is  the  flexural  rigidity,  C = Eh, 
where  E is  Young's  modulus,  and  v is  Poisson's  ratio,  then  the  basic 


differential  equations  are 

(4a)  (D/ a)[  3"  + 3'  cot  Z - (cot2§  + v)  3]  = -'ll  - -|pa“3  + \jf3  cot  £ , 

(4b)  ij<"  + i|r’  cot  Z -(cot2§  - v)  t = Ca3  - |CaP2cot  Z , 


where  the  prime  signifies  differentiation  with  respect  to  § , Equations 

(4a)  and  (4b)  are  deduced  from  [12,  Eqs.  (66)  and  (67)]  by  setting 

Pjj  = Py  = 0 and  neglecting  all  terms  on  the  left  which  do  not  appear 

in  the  classical  Reissner-Meissner  equations  and  all  those  on  the  right 

which  do  not  have  corresponding  terms  in  the  shallow-shell  approximations 

to  these  equations  [12,  Eqs.  (72)  and  (73)]. 

Equations  (4a)  and (4b)  already  permit  the  determination  of  P . 

L 

However,  it  is  desirable  to  be  able  to  compute  the  buckled  shape  of  the 

shell,  from  3 and  Ji  , For  that  purpose  we  consider 

Eqs.  (63),  (68),  and  (6y)  of  [12]  and  obtain  after  setting  P = P = 0, 

H V 

the  dimensionless  displacements 


(5a) 


(1-v) 


Pa  . P j/’ 

2 Eh  Sln  ^ + Eha 


sin 


-i  - V — - — COS 

Eha 


Pa  q 2 - r 

rr-  & sm  Z cos  Z~ 


_1_ 


v P sin  ^ 
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t 


(5b)  w(f) 


^ £ j' ...  s « - j‘ a 2 Eh 


(i-v)  £p_)(3  cos  £ d£ 


§ 


- v J\  IL  sin  5d5+  L imCOR  ?d? 


v J 1^-0  sin“§  cos  § d§  “ i J p2sin  5 

^0  ^0 


§ 5 

+ J ^ 3 sin  § d§  - J_  Ihl  3 COt  5 COS  5 d? 

T)  % 


5 t 5 

+ v r 4—  3 COS  § d§  + V J 1^-  p2sin  § cos2§  d§  , 


i”.  Eh  a 
T) 


where  0 < ^ ^ n is  that  value  of  §,  to  bo  chosen  ad  libitum,  for 
which  the  axial  displacement  vanishes,  i.e.,  W(§Q)  = 0.  The  first  terms 
on  the  right  of  (5a)-(5b),  respectively,  are  the  respective  components  of 
the  membrane  contraction  with  the  normalization  indicated. 

It  should  be  noted  that  all  cubic  terms  in  Reissner's  formulas  have 
been  deleted  so  that  this  is  a "quadratic"  theory  (see  Section  7 below 
for  further  comment  on  this  point).  The  modifications  of  the  equations 
necessary  to  study  imperfections  are  indicated  in  Section  7. 


i 

i 


i 


1 

* 
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3 . Method  of  Solution  of  the  Equations 

The  key  to  the  solution  of  (4a) -(4b)  and  the  subsequent  evaluation 
of  (5)  is  the  observation  that,  the  common  differential  operator 


(6)  ()"+()'  cot  S - ( ) cot2S 

on  the  left  of  (4a) -(4b)  can  be  written 

()"+()'  cot  § - ( )(csc2§  - 1)  = L(  ) + 1(  ), 

where 

L(P1(cos  p)  = -nOi+DP1  (cos  p, 
n n 

Pn(x)  being  the  associated  Legendre  function  of  order  one  and  degree  n 

[15,  Chapter  XV,  §15.5].  Thus  P1 (cos  p is  an  everywhere  finite  eigen- 

n 

function  for  (6)  with  eigenvalue  1 - n(n+l) . 

Two  facts  about  the  Legendre  polynomials  P (x)  and  associated 

n 

Legendre  functions  of  first  order  should  be  noted.  One  has,  with 
x = cos  p 

(7a)  P1(x)  = (1  - x2)^  dP  (x)/dx  , 

n n 

(7b)  (d/dpp  (cos  P = -(sin  PdP  (cos  p/dx  = -P^(cos  p, 

n n n 

(7c)  [d2Pn(cos  p/d^2]  + (cot  p(d/dpP^  (cos  p = -n(n+l)Pn(cos  p, 

(7d)  (1  - x2)[d2P1(x)/dx2]  - 2x(d/dx)P1(x)  - [1/(1  - x2)]P1(x) 

n n n 

= -n(n+l)P1(x)  . 

n 1 

For  technical  reasons,  which  will  become  apparent  later,  it  is  con- 
venient to  introduce  x = cos  § as  independent  variable  in  (4a) -(4b)  to 
obtain  [ it  should  be  understood  that  we  are  writing  3(x)  = fKcos  p 
rather  than  3(P,  etc.] 


(8a) 


5[(1  - X2)  sf|  - 2x  - (— + vVl  = —■  B + HP  . 

aL  dx2  dx  V 1-X2  AJ  2 (1-xV 

2 


(8b)  (1  - x2)  O - 

dx 


„ 2x  U 
2 dx 


Mr  - v>  - 


1— x 


Ca  fl2  x 

CaP  “ P o i * 

2 (1-xV 


\ 


sinco  ( 6 ) becomes 
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(9)  (1  ~ x2)  (d2/dx2>  ( ) - 2x(d/dx)(  ) - [x2/(l  - x2)](  ). 

Since  Pn(x) , n = 1,..,,  form  a complete  orthogonal  sot  on  (-1,1), 
it  seems  reasonable  to  use  them  to  find  approximate  solutions  of  (8a) -(8b) 
for  prescribed  p by  the  spectral  method.  We  set 

s i 

(10a)  8 = £ A P (x) 

. n n 
n=l 


(10b) 


♦ * £ B P (x) 

, n n 
n=l 


for  fixed  s,  substitute  in  (8aM8b) , oxpand  the  right  sides  out  in 
P*(x),  n = 1,...,  retaining  only  the  first  s terms,  and  then  compare 
sides. 

As  a first  trivial  but  vital  application  of  the  method  we  determine 

the  classical  linear  buckling  load  p . Set  8 = AP1  and  & = BP1 

crxt  n T n 

in  the  equations  obtained  from  (8a)-(8b)  by  ignoring  nonlinear  terms, 
and  obtain 

(D/a)[l  - n(n+l)  - v]  = -B  - §pa2A,  [l  - n(n+l)  + v]B  = CaA  . 

On  eliminating  B and  assuming  A ^ 0,  one  obtains 
3 


(ID 


2E(h/a)~  , / , 2Eh  1 

P = — - — [ n(n+l)  - 1 + v]  + 


a n(n+l)  - 1 - v 


12(1- v ) 

Differentiating  (ll)  with  respect  to  n and  equating  the  result 
to  zero  give 

- v2i 


(12) 


n(n+l)  - 1 - v = 2[3(1  - v )]2(a/h)  . 


Substituting  in  (11)  yields 


(13)  p 


2E 


crit 


[3(1-V)]2 


(if 

2.  , o'  \a  / 


il  + 


2[  3(1-  \?)  ] s a' 


1 £) 


2E 

[ 3 (1 -v2) ] 


I C-)  <-)  ■ 


I 


-f)~ 


where  L(h/a)  is  tho  parenthesis. 
Returning  to  (10) , wo  define 


(14a) 


(-t.m.n)  = J 


■7)  T-  P^(x)P1(x)p\x^dx 
v m n 


-1  (1-x  ) 

1 „ dP.  dP  dP 

r 2.  t ra  n 

J dx  dx  dx 


(14b) 


j = J (P  ) dx  = 
n J , n 


1.2  2n(n+l) 


-1 


2n+l 


and  carry  out  the  spectral  method,  first,  by 
substituting  (10)  in  (8b)  and  comparing  sides  to  obtain 

s 

(A)  B = -fCa/[  n(n+l)  -1-  lA  + [Ca/2[  n(n+l) -1-v]  } E [ (■t,m,n)/j  ]A  A , 

n n # . n t,  m 

•t,  m=l 

n = l,...,s,  and  then  by  substituting  (10)  in  (8a)  to  obtain 
„ s 


(B) 


( gpa  - (D/a)[n(n+l)  - l+v])A  = -B  + E [ (£,m,n)/j  ] A B , 

n n i*=l  n m 


n — 1 , • . . , s * 


We  introduce  the  dimensionless  variables 


(15) 


* , 2 
B = B /p  a 
n n ent 


p’  - p/Pcrit  ' P”2/Pcrlt  »2  • 


2 

Dividing  equations  (A)  and  (B)  by  p .a  , multiplying  (B)  by  two,  and 

crit 

using  (15)  give 

[3(1-v2)1^ 

2L(h/a)[  n(i 

<A,)  [3(1- v2)]* 

+ ..-I.  r 


B*  _ L 3(l-v~)  J c a 

n 2L(h/a)[  n(n+l) -1- vj  h n 

s 


[ 3(1-  V ) 1 * a J r(^,m,n)-[, 

4L(h/a)[  n(n+li-l-  vj  h L j i*  vi 

in—  X n 


jp*  - H } A 

1 a 4[3(l-v2)]%(h/a)  j n 

~(-t,  m,n) 


(B')  _ s 


+ 2B*  - 2 E r1-^-  1 AV  = 0 . " = 1. 

" t,u*l  L Jn  J 1 m 


1 


To  compute  the  shape  o.C  tho  shell  corresponding  to  any  solution 
* + * 

P , , I5^,.,.,Il  of  (A')  and  (!>')  wo  represent  tlio  middle 

surface  in  dimensionless  form  as 

(1C)  r 10  u _ u 

— = — + — = sm  q + — 

a a a a 

z 

z 0 w _ w 

— = — + — = - cos  | + — 

a a a a 

and  write  in  Eqs.  (5a)  and  (5b), 

(17)  pa  h 2L(h/a)p* 

Eh  “ V r-- 2. 

V§(l-v  ) 

* ,h,  2L(h/ a) 

Eh  a a 5 — * 

y£(l-v  ) 

With  the  last  definition  one  has 

**  =,  IB*  P1(cos  5)  , 

, n n 
n=l 

and  if  one  substitutes  this  and  the  truncated  series  for  B in  (16) 
and  (17)  then  numerical  integration  using  the  trapezoidal  rule  yields 
the  shape.  Tho  shape  computation  was  effected  by  a simple  auxiliary 
FORTRAN  program.  The  value  of  was  taken  to  be  n,  i.e.,  the  north 

pole  was  taken  to  be  fixed.  The  advantage  of  such  a normalization  will 
be  apparent  in  a moment. 

In  tho  program  SPHERE  (see  Section  4)  used  to  solve  (A')  and  (B')  it 

was  necessary  to  incorporate  a subroutine  which  for  each  solution 

* 

P , A^ , . . . computes  some  convenient  deflection  parameter.  At  first 
glance  (and  in  [ll])  it  seems  reasonable  to  use  the  south  polar  deflection 
under  the  assumption  that  the  equator  is  fixed  ( ^ = tt/2).  In  hindsight, 
if  the  above  mentioned  auxiliary  program  had  been  used  as  the  subroutine, 
as  is  possible,  then  the  fixed  equator  hypothesis  would  have  been  as 


-li- 


en sy  to  implement.  However,  the  original  polar  deflection  subroutine 

jj;  ]|( 

uses  the  coefficients  A^,...,A  ; B^.,.,1^  directly,  and  there  as  will 
bo  seen  in  Section  C,  tho  simplification  brought  about  by  setting  ^ = n 
is  enormous.  We  then  docided  that  half  the  resulting  deflection  is  the 
useful  parameter  for  tho  symmetric  deformations  (see  Section  5) . Since 
the  resulting  formulas  are  in  SPHERE  and  may  be  of  independent  interest, 
wo  reproduco  them  here. 

We  define 


(18) 


A(m,n)  = A(n,m) 


L 

2m+l 


m ^ n 

> 

m = n 


p(n) 


q 
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2L(h/n) 

-^(1-v2) 


n = 0 (2) 
n = 1 (2) 

asymmetric  deformation 
symmetric  deformation. 


Thon  the  polar  deflection  parameter  v is  given  by 
(19)  qv  = 2^21  = (l-v)Kp*  + (l-(— Kp*)  E (l-(-l)n+1)A 

A . Il 
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This  formula  will  be  established  in  Section  6.  In  practice  only  the 
linear  and  the  first  non-linear  terms  make  significant  contributions. 


4 . Numerical  Methods 

For  given  a/h,  v,  and  s numerical  solution  of  the  equations  (A') 
and  (B’)  was  performed  by  means  of  the  program  SPHERE,  written  for  this 
purpose  but  adaptable  to  other  sets  of  non-linear  equation.  SPHERE  is 
modeled  on  the  algorithm  proposed  by  E.  Polak  [ 10] . The  essence  of  the 
algorithm  is  a combination  of  the  secant  method  with  the  method  of  local 
variations.  We  found  it  to  be  very  powerful. 

Already  in  our  first  version  of  SPHERE,  a PL/I  double  precision 


program,  we  omitted  certain  convergence  accelerating  features  at  the  end 
of  Polak 's  algorithm.  It  should  bo  pointed  out  that  there  was  little  hope 
that  Polak 's  hypothesis  on  invertibility  of  the  Jacobian  would  bo  verified 
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i 
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In  our  case  and,  thus,  little  hope  of  global  convergence.  Indeed,  in 
continuation  of  solutions  near  the  first  relative  minima  of  the  pressure- 
deflection  graphs  we  Often  encountered  suggestions  of  possible  zero 
Jacobians  in  the  form  of  loss  of  significance  messages  from  the  linear 
inversion  subroutine.  Correspondingly,  we  had  little  success  in  using 
SPHERE  to  find  many-mode  solutions  starting  from  wild  guesses.  Accord- 
ingly the  necessity  for  good  starting  and  continuation  procedures 
manifested  itself  early.  We  shall  describe  some  below. 

Equations  (A’)  and  (B')  were  cast  in  the  functional  form  G(A)  = 0, 
where  G(A)  is  a vector  with  s components,  by  putting  all  terms  on  the 

left  side  in  each  equation  of  (B1)  and  then,  for  each  attempted  solution 

* * , t. 
vector  A^,...,As,  substituting  the  vector  B^,...,Bs  computed  from  (A  ) 

into  (B').  As  a criterion  for  a "solution"  we  used  |(g (A)  ||  £ e , where 

the  norm  is  the  usual  vector  norm,  and  e was  a predetermined  small 

number,  usually  10  ® . 

The  coefficients  (-t,m,n)  required  in  G(A)  = 0 for  a given  s 
were  computed  by  a program  TABLE  on  the  basis  of  equation  (27)  in  Section 
6 and  put  on  a disk  which  was  read  into  core  when  SPHERE  was  run.  Because 
of  the  large  size  of  the  three-way  array  (£,m,n)  for  realistic  s - 
despito  the  symmetry  and  vanishing  properties  of  Section  6,  which  reduce 
64,000  to  3270  for  s = 40,  for  example  - and  the  initial  demand  for 
double  precision  a linearization  routine  for  both  TABLE  and  SPHERE  was 
written  (by  James  Korenthal).  However,  it  developed  that  single  precision 
was  quite  adequate  and  that  the  linearization  was  quite  time  consuming 
so  that  the  (-t-,m,n)  could  be  road  in  directly  with  tolerable  core 
requirements  and  very  substantial  savings.  Ultimately  a still  simpler 


and  faster  FORTRAN  version  of  SPHERE  was  prepared.  Copies  of  all  programs 
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dcscribed  in  this  paper  are  available  by  request  fi\>in  the  authors. 

All  <ii sou. ss ion  oi'  accuracy  and  start  lng  and  continuation  procedures 
hinges  on  the  notions  of  critical  motto  and  supercritical  mode.  Here, 

a mode  is  one  the  A's,  the  n-th  mode  is  A , and  n is  the  index  of 

n 

that  mode.  Then  a mode  is  critical  for  a Riven  value  of  a/h  if  its 

index  minimizes  the  expression  (11)  for  p obtained  from  the  assumption 

* 

of  linear  buckling  or,  equivalently,  the  expression  for  p obtained  by 
dividing  (11)  by  (13),  where  a/h  is  the  given  value.  The  mode  whose 

index  differs  from  the  minimizing  index  by  one  and  whose  associated  value 

* 

of  p is  the  next  lowest  is  also  called  critical.  Thus  to  each  value 
of  a/h  is  associated  uniquely  a pair  of  critical  inodes  whose  indices 
are  consecutive  integers,  in  particular,  one  odd  and  one  even.  For  the 
two  values  of  a/h,  100  and  200,  handled  explicitly  in  this  paper  one  has, 
for  a/h  = 100,  the  values  p*  = 1.019939,  1,003290,  1.000465,  1.009291 
associated  with  the  indices  16,  17,  18,  19,  respectively,  all  other 


indices  giving  higher  values,  so  that  the  seventeenth  and  eighteenth  modes 

* 

are  critical  with  the  even  mode  having  lowest  p . Similarly  for 
a/h  = 200  one  has  1.004888,  1.000174,  1.001709  corresponding  to  24, 

25,  26,  respectively,  so  that  modes  25  and  26  are  critical  with  the 

* 

odd  mode  having  lowest  p . A supercritical  mode  for  given  a/h  is  one 
* whose  index  is  greater  than  those  of  the  critical  modes  for  that  a/h. 

The  critical  modes  play  a role,  first,  in  starting  procedures.  The 

* 

basic  starting  strategy  to  find  a solution  of  G(A)  = 0 for  given  p 
and  a/h  is  to  solve  the  smallest  number  of  equations  possible,  i.e., 
uso  the  smallest  number  of  modes  possible,  and  then  extend  this  solution 
to  one  with  more  modos  by  using  its  entries  as  the  initial  entries  of  a 
guess  for  a higher  number  of  modos,  witli  the  remainder  of  the  entries 
being  zeros.  Wo  call  this  strategy  extending  a solution  by  zeros.  It 

Li  - l , „ „ J 
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is  hampered,  howovez-,  by  the  pccular  fact,  discovered  in  the  course  of 
tho  numerical  work,  that  it  Uoos  not  work  unless  the  number  of  modes  or, 
equivalently,  tho  index  of  tho  last  (highest  indox)  modo  of  tho  solution 
to  be  extended  is  greater  than  tho  indices  of  tho  critical  modos  for 
that  a/h.  Moro  precisely,  only  solutions  already  containing  super- 
critical modes  can  be  extended  by  zeros  to  obtain  solutions  with  (even 
more)  supercritical  modes. 

"Subcritical"  solutions  can  often  be  extended  by  zeros  to  "longer 
subcritical"  solutions,  but  these  cannot  be  extended  further  and  represent 
spurious  approximations  to  solutions,  analogous  to  the  extraneous  solu- 
tions to  the  difference  equations  sometimes  obtained  when  applying 
finite-difference  methods  to  differential  equations. 

Accordingly  it  is  necessary  to  start  by  assuming  small,  quite  un- 
physical, values  of  a/h  such  as  .5  for  which  A is  supercritical. 

* 

For  fixed  p (for  example,  a little  less  than  1,  if  one  is  starting  a 

pressui-e-def  lection  graph  of  the  postbuckling  regime)  one  solves  for 

A , A by  starting  with  any  reasonable  initial  vector,  say,  (1,0)  - 
1 / 

an  even  more  surefire  technique  is  given  below.  The  resulting  solution 
is  then  extended  by  zeros  with  a sufficiently  large  number  of  terms  so 
that  the  last  index  is  supercritical  for  a higher  value  of  a/h,  say,  5. 
The  resulting  solution  is  then  continued  by  letting  a/h  vary  from  the 
starting  value  to  the  new  value.  Tho  new  solution  is  then  extended  by 
zeros  until  supercritical  for  a higher  value  of  a/h,  and  so  on  until  a 
solution  at  the  desired  a/h  is  obtained.  This  is  called  zig-zagging. 

It  proved  to  be  more  efficient  to  zig-zag  in  reasonably  small  increments 
of  a/h  rather  than  to  attempt  to,  say,  extend  the  initial  solution  to 
be  supercritical  for  tho  final  valuo  of  a/h  and  then  to  continue  to 
that  valuo.  Also,  it  was  found  to  bo  much  easier  to  oxtend  solutions  by 
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zcros  or  continue  solutions  in  the  parameter  a/h  for  p noar  1 

* * 

and  small  deflections  rather  than  for  p near  p and  the  corresponding 

Ij 

lai'i'o  deflections. 

A pressure-deflection  graph  is  not  monotone  in  principle,  and  the 
approximations  to  it  obtained  by  truncating  the  expansions  at  too  small 
an  index  are  often  even  more  oscillatory.  Hence,  the  vital  feature  of 
continuation  by  changing  one  of  the  A's,  typically  A or  A for 

JL  £. 

asymmetric  or  symmetric  solutions  (sec  Section  5),  respectively,  and 
* 

solving  for  p as  one  of  the  unknowns  was  built  into  the  program.  Once 
this  was  done  it  was  found  that  it  was  much  more  efficient  to  use  this 

* 

method  of  continuation  even  on  monotone  branches  so  that  solution  with  p 

♦ 

fixed  and  continuation  by  varying  p wore  used  ultimately  only  in  the 
starting  procedure  mentioned  above.  It  was  found  that,  in  general,  al- 
though not  always,  the  polar  deflection  is  a monotone  increasing  function 
of  A or  A . 

X 

For  very’  large  values  (say,  > 15  radii)  of  the  polar  deflection  and 
correspondingly  large  A's  it  is  necessary  to  relax  the  demand  on  the  norm 
of  the  residuals  by  several  orders  of  magnitude  in  order  to  continue 
computing  the  graph  efficiently  since  the  residuals  are  extremely  sen- 
sitive to  small  changes  in  the  A's  when  the  first  few  of  the  latter  are 
large,  as  one  sees  by  inspecting  (A')  and  (B').  However,  any  one  solution 
can  be  refined  and  the  effect  on  tho  graph  is  visually  imperceptible. 

It  should  be  noted  that  if  one  distinguishes  symmetric  from  asymmetric 
solutions  to  (4)  (see  Section  5)  then  the  associated  truncated  expansions 
can  bo  identified  easily.  In  fact,  as  noted  there,  tho  coefficients  A 
and  Li^  witli  odd  n vanish  for  a symmetric  solution.  This  rcducos  the 
acl ual  number  of  variables  for  even  s to  s/2.  SPHERE  is  programmed 


to  tako  advantage  of  this  by  handling  only  variables  with  even  indices 
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upon  receipt  of  the  proper  input,  effecting  a very  substantial  saving 
in  computing  time. 

Finally  wo  outlino  a foolproof  procedure  for  solving  G(A)  = 0 
when  s = 2.  Wo  note  that  by  properties  of  the  (-t,m,n)  obtained  in 
Section  6,  (1,1,1)  = (2,2,1)  = (2,1,2)  = (1,2,2)  = 0 and  (2,1,1)  = 
(1,2,1)  = (1,1,2).  For  s = 2 (A')  and  (B')  become 


B?  = + |(2-,1,1)A2A1) 


1 2L(h?a)  (1-vY 

T 


B2  - aiLrili-a  d)(-A2  + §4<1-1-2)ai  + h(2’2-2)Al) 


0>*  - (t) 


(1  + v) 


v6(l-v2)L,(h/a) 


) A1  = _2B1  + 


2 + 


0>*  - C-) 1 


4, 

<1  + v) 


r 


)a  = -2B*  + |(1,1,2)A1B*  + |(2,2,2)A0B*  • 


2 2 


4/^6  (1-v  )L(h/ a) 


If  one  sets  A^  = 0,  B^  = 0 so  that  one  is  seeking  a symmetric  solution, 
then  one  sees  that  the  first  and  third  equations  are  identically  satisfied. 
One  could  then  substitute  the  second  in  the  fourth  to  obtain  an  explicit 
cubic  for  A which  could  be  solved  numerically,  but  that  is  not  necessary 
if  one  is  using  SPHERE  since  the  use  of  initial  guess  vectors  with  the 
first  component  zero  is  equivalent  to  solving  that  cubic  by  a secant- 
local-variation  method.  If,  however,  one  desires  an  asymmetric  solution, 
so  that  A^  ^ 0,  then  one  finds  that  substitution  of  the  first  two 
equations  in  the  third  leads  to  a cubic  equation  in  A and  A every 
term  of  which  contains  A^  to  the  first  or  third  power.  Cancellation 


of  A^  leads  to  a quadratic  which  can  bo  solved  explicitly  for  as 


follows . 


32  (5-v) 
5 (1-v) 


_1 

(1,1,2)  2 


fX(h/n)  (1-v) 
i 2 
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Again  this  equation  could  be  substituted  along  with  the  first  and  second 
equations  into  the  fourth  equation  to  obtain  an  explicit  cubic  in  A^ 
whose  solution  together  with  this  equation  would  solve  G(A)  = 0 for 
s = 2.  Hut,  again,  it  is  more  suitable,  if  one  has  guessed  A^,  to 
calculate  A^  from  the  preceding  equation  and  use  the  resulting  A^,A^ 
as  input  for  the  program. 

It  should  be  mentioned  that  the  discussion  of  imperfection- 

sensitivity  at  the  end  of  Section  7 suggests  another  starting  procedure. 

What  makes  starting  so  difficult  for  (A1)  and  (B1)  is  that  they  are 

homogeneous  so  that  they  possess  a solution  whose  entries  are  all  zero. 

By  inserting  a non-zero  imperfection  parameter  e in  the  equations  at 

n 

the  end  of  Section  7 one  obtains  inhomogeneous  equations  of  which  non- 

* 

trivial  solutions  which  continue  the  trivial  solution  for  p = 0 are 

♦ 

easily  found  for  small  p . Continuation  to  the  relative  maximum 

(or  maxima)  and  then  to  the  downward  sloping  part  of  the  graph  followed 

by  continuation  in  e as  e -•  0 would  yield  non-trivial  solutions  of 

n n 

the  equations  for  the  perfect  shell  which  could  then  be  continued  for 
smaller  or  larger  deflection. 


5.  Numerical  Results 

As  a preliminary  we  note  that  replacing  5 by  tt  - § in  (4)  while 
replacing  3,ijf  by  -3,  -t»  respectively,  shows  that  if  3(D,  l|l(0  are 
a solution  then  -3(TT_§)>  -V(TT-§)  are,  too.  Equivalently,  (with  tho 


same  abuse  of  functional  notation  remarked  on  in  Section  2) 


-3(-x),  -'K-x)  arc  also  a solution  of  (8).  From  the  remarks  at  the 
beginning  of  Section  2 it  is  clear  that  this  means  the  physically 
obvious  fact  that  if  one  has  a buckled  shell  then  the  mirror  imago  in 
the  equatorial  plane  of  the  unbuckled  shell  is  also  a buckled  shell. 

In  particular,  a solution  0,  i)i  for  which  0(§)  = -0(Tr-§)  , ^(§)  = -<ji(n-5) 
or,  equivalently,  3(x)  = -0(-x),  (J'(x)  = -^f(-x)  will  be  called  symmetric 
because  it  manifestly  represents  a buckled  shape  symmetric  with  respect 
to  the  equatorial  plane.  By  the  preceding  remarks  the  asymmetric 
solutions  come  in  pairs,  one  being  the  mirror  image  (in  the  equatorial 
plane)  of  the  other.  In  terms  of  the  truncated  expansions  of  0 and  (jf 

1 2 h 

one  observes  that  P (x)  = (1-x  ) EdP  /dx  and  the  fact  that  P (x)  is  even 

n n n 

or  odd  according  as  n is,  together  with  0(x)  = -0(-x),  ^i(x)  = -t(-x) 

* 

imply  A = B = 0 for  odd  n,  that  is,  a solution  is  symmetric  if  and 
n n 

only  if  terms  with  even  indices  only  appear  in  the  expansions.  For  this 
reason  symmetric  solutions  are  referred  to  as  even  in  our  programs  al- 
though, in  fact,  they  are  odd  as  functions  of  x.  Similar  reasoning 
shows  that  one  of  the  paired,  mirror-image, asymmetric  solutions  differs 
precisely  from  the  other  in  the  signs  of  the  terms  of  odd  index.  The 
presence  of  terms  with  odd  indices  in  the  expansions  of  asymmetric 
solutions  led  us  to  call  them  odd -even  in  the  programs  although  they 
are  neither  odd  nor  even. 

For  all  the  values  of  a/h  used  in  our  work  we  found  for  any  even 

supercritical  s three  solution  branches  of  G(A)  = 0 issuing  from  the 
% 

vicinity  of  p = 1,  A^  = ...  = A^  = 0,  i.e.,  the  (unbuckled)  membrane 

solution  at  the  classical  linear  buckling  pressure.  By  a solution 

* 

branch  wo  moan  a one-parameter  family  of  values  of  p and  associated 
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non-trivinl  vectors  A each  of  which  satisfies  G(A)  - 0 for  t ho 

* * 

associated  value  oi'  p . The  parameter  is  either  p itself  or  one 

of  the  A's,  usually  A or  A . Ono  of  the  branches  consists  of  sym- 

1 2 

* 

metric  solutions  and  issues  from  the  point  for  which  p is  the  value 
associated  with  the  critical  mode  of  even  index  for  the  given  a/h,  and 
the  other  two  issue  from  the  corresponding  point  for  the  critical  mode 

of  odd  index,  any  one  of  the  two  representing  the  mirror  images  of  the 

* 

solutions  of  the  other.  When  the  value  of  p and  the  value  of  the  polar 
deflection,  which,  wo  recall,  is  the  right  side  of  (18)  for  asymmetric 
solutions  and  half  that  for  symmetric  ones,  are  plotted  as  ordinate  and 
abscissa,  respectively,  one  obtains  two  pressure  deflection  graphs,  one 
for  the  symmetric  solution  branch  and  one  for  the  two  mirror  image  asym- 
metric solution  branches  since,  as  ono  verifies  from  (5b),  the  pola- 

def lection  is  the  same  for  mirror  images.  We  ignore  the  prebuckling  linear 

* 

membrane  pieces  of  the  graphs  since  they  play  no  role  in  determining  p . 
However,  if  imperfection-sensitivity  were  studied  in  the  way  sketched  in 
Section  7,  then  the  corresponding  parts  of  the  graphs  would  be  quite  im- 
portant and  would  require  much  larger  scales  than  those  of  Figures  2 and  3 
here . 

We  take  the  lessor  of  the  absolute  minima  of  the  two  graphs  in  question 
* 

for  given  a/h  to  be  p . Now  for  a/h  = 100,  200,  the  two  cases 
handled  in  extenso,  the  two  graphs,  symmetric  and  asymmetric  solutions, 
in  each  case  are,  most  remarkably,  virtually  coincident  and  certainly  so 
when  graphed  on  any  reasonable  scale.  For  that  reason  Figure  2 actually 


is  two  graphs,  and  the  same  is  true  of  Figure  3.  The  second  interesting 
fact  is  that  for  botli  100  and  200  the  asymmetric  minimum  proved  greater 
than  the  symmetric  minimum  but  only  by  at  most  ono  digit  in  tho  last 
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significant  figure  l’otaincd.  In  [l]  for  a/h  = 9.13  (their  Figure  lb) 

it  is  also  true  that  the  symmetric  minimum  is  less  than  the  asymmetric 

minimum.  On  the  hypothesis  that  this  is  always  true  for  9.13  £ a/h  £ 200 

wc  have  included  two  other  sets  of  data  in  Tablo  1,  a table  of  lower 

critical  pressures.  One  is  a run  for  a/h  - 21.4  we  made  for  symmetric 

solutions  only  and  tho  other  for  a/h  = 91.3  is  takon  from  Figure  9 of 

[ 1]  by  graphical  estimation.  Other  data  included  in  Table  1 for  tho 

sake  of  comparison  are  the  values  of  k and  P , the  parameters  used 

L 

in  [ 1 ],  corresponding  to  our  a/h  and  p*  , respectively, as  well  as  the 

L 

deflection  parameter  A of  [l]  for  the  relevant  solutions  and  our  polar 
deflection  parameter  v for  our  solutions  and  those  of  [ l]  where  graph- 
ical  estimation  permitted  us  to  determine  it.  A is  a root  mean  square 
deflection,  which  wo  did  not  compute  for  our  solutions. 

We  call  attention  to  another  unanticipated  result  of  our  computations. 
Figure  4 is  the  shape  of  the  asymmetrically  deformed  shell  for  p = p*  at 

li 

a/h  = 100,  but  it  is  also  the  shape  of  the  corresponding  shell  for 
a/h  = 200,  within  the  tolerance  of  the  graphical  realization.  There  is 
a hint  here  of  an  asymptotic  shape  as  a/h  -*  ® . if  one  regards  the  shell 
as  an  inoxtensible  membrane  then  a classical  theorem  of  differential 
geometry  in  the  large  implies  that  the  shell  is  rigid,  i.e.,  admits  no 
(differentiable)  bending.  The  only  possible  deflected  shape  is  obtained 
in  the  obvious  way  by  violating  differentiability  and  sawing  off  a cap 
and  reversing  it.  Now  a careful  inspection  of  Figure  4 discloses  that 
tho  indentation  resembles  a piece  of  sphere  of  radius  a subject  to 
increasingly  strong  bending.  But  in  the  absence  of  elasticity  theory 
there  is  no  a priori  reason  to  select  one  sizo  of  cap  over  another,  where- 
as Figure  4 shows  that  the  non-linear  elasticity  theory  embodied  in 


Re issuer's  equations  predicts  a definite  size,  a central  angle  of  a 
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liLtlo  less  than  120  for  a/h  = 100,  200  and  possibly  asymptotically. 

Two  final  points  of  discussion:  accuracy  and  physical  significance. 

The  figures  in  Tahiti  1 Lor  a/h  100  and  200  arc  accurate  at  least  to 
the  number  of  significant  figures  shown.  This  was  determined  by  in- 
creasing the  number  of  modes  until  the  figures  stabilized.  The  results 
shown  used  40  modes  for  a/h  = 100  and  60  for  a/h  = 200.  The  last 
probably  approaches  the  limit  of  the  present  codes  within  practicable 
computer  times.  More  modes  and  greater  a/h  require  improved  codes 

(see  Section  7).  As  a sample  of  our  results  Table  2 shows  A , . . . jA^q  ; 

* * 

B , ,..,B  (rounded  to  four  figures' to  save  space) for  a/h  = 100, 

♦ 

p = .0674,  which  is  the  asymmetric  minimum. 

As  far  as  physical  significance  is  concerned  the  likelihood  or  lack 
of  it  of  realizing  deflections  of  the  magnitude  shown  in  Figure  4 is 
discussed  in  Section  7.  However,  the  physical  usefulness  of  the  values 

jJ* 

of  P in  Table  1 for  the  thin  shells  (a/h  S 91.3)  should  not  be  dis- 

L "■ 

missed  lightly.  Discussions  of  imperfection-sensitivity  indicate  that 

* 

the  snap-through  load  can  be  degraded  substantially  toward  p^  . 


. — 


** 

-3 

15.7  % 

3.4  $ 

9.13 

10  J 

.0105 

2.5 

21.4 

1.82X10-4 

12.4 

-3  4 

3.52X10  T 

3.91  ^ 

- 

** 

91.3 

io-5 

1§ 

7.2  %TS 

-44  § 

4.8X10 

- 

1.2 

-6 

6.73  % 

-4 

100 

8.3>d0 

4.08  XL0 

1.41 

“ 

200 

2.08  X10“6 

4.9  % 

1.5XL0-4 

1.4 

— 

Results  from  [ 1]  . 


'Computed  from  symmetric  solutions  only. 
^Estimated  from  graphs  in  [l]. 
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TAIiLE  2 


EXPANSION  COEFFICIENTS,  a/h  = 100,  v = 1.41,  p*  = .0674.  p*  = 
" L 


A*  s (read  by  rows;  3.063  - 1 = 3.063  X 10  , etc.). 


3.063-1 

3,706-1 

2.982-1 

1.400-1 

-1.858-2 

-1.004-1 

-9.791-2 

-2.906-2 

4.070-2 

6.475-2 

3.749-2 

-1.000-3 

-3.911-2 

-3.306-2 

-4.589-3 

2.011-2 

2.398-2 

9.561-3 

-7.882-3 

-1.494-2 

-9.381-3 

1.319-3 

8.077-3 

7.192-3 

1.462-3 

-3.669-3 

-4.740-3 

-2.156-3 

1.208-3 

2.753-3 

1.912-3 

-3.514-5 

-1.396-3 

-1.380-3 

-3.949-4 

5.859-4 

8.719-4 

4.585-4 

-1.693-4 

-5.175-4 

-1 .462-2 

-1.201-2 

-6.366-4 

1.110-2 

1.474-2 

8.111-3 

-3.290-3 

-1.081-2 

-9.419-3 

-1.305-3 

6.499-3 

8.153-3 

3.446-3 

-2.888-3 

-5.743-3 

-3.724-3 

5.447-4 

3.353-3 

3.002-3 

5.751-4 

-1.605-3 

-2.030-3 

-8.863-4 

5.649-4 

1.200-3 

8.022-4 

-4.930-5 

-6.201-4 

-5.882-4 

-1.505-4 

2.667-4 

3.756-4 

1.886-4 

-7.707-5 

-2.125-4 

-1.601-4 

-1.033-5 

1.054-4 

1.173-4 

4.803-5 
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6.  Legendre  functions  and  the  coefficients  (-1,01, n) 


Wo  need  the  recursion  formulas  [15,  §15.  21] 


xdP  /dx  = nP  + dP  /dx 
n n n-1 


(1-x  )dP  /dx  = nP  , - nxP 
n n-1  n 


(22)  xP  = [ (n+l)/(2n+l)]P  + [n/(2n+l)]P  , 

n n n-1 


(23)  dP  /dx  = dP  ./dx  + (2n-l)P  , . 

n n-2  n-1 

From  (20)  and  (21)  we  deduce 

(24)  (l-x2)dP  /dx  = [n(n+l)/(2n+l)] (P  -P  ,) 

n n-1  n+1 

From  (22)  by  induction  we  deduce 

(n+p (n))/2 

(25)  dP/d„  . Z <“-l-1-2P<"»P2J-1-p(n) 

j * 


Finally  from  (20)  and  (25)  we  deduce 

((n-1) +p (n-1 1/2 

(26)  xdP  /dx  = nP  + £ (4j-l-2p(n-l))P 

n j=l 


2j-l-p(n-l) 


(n -p  (n))/2 

= nPm  + £ (4j-3+2p(n))P 


2j-2+p(n)  ’ 


since  p (n-l)=  1-p (n) . 

Now  apply  (25)  with  n = t,  (24)  with  n = m,  and  (26)  to  obtain 

1 o dpJ  dp  dp 

(t,n,n)  a J x(l-x  ) ——  - — — dx 

dx  dx  dx 

_ m(m-t-l)  (,tlp(<t)y2(4i-l-2p(-0)fnrP(2i-l-p(<t).m-l.n) 


(n-p(nl/2 

- P(2i-l-p(£)  ,m+l,n)]  + £ (4,i-3+2p(n)) 

j=l 

X [P(2i-l-p(/.)  ,m-l,2J-2+p(n))  - P(2i-l-p(i)  ,nn-l , 2j-2+p(n)  )]  } , 
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JL 

where  P(-t,m,n)  = f P-P  P dx  . 

" . -t  m n 


Now  according  to  a result  of  Ferrers  and  Adams  [15,  p.  3311, [13,  p.306] 


(28) 


JL 

f P.P  P dx  = 
0 L m n 


_2 A_£  e— t)  A ( G-m)  A ( e-n) 

2Gil 


-1 


A ( g) 

if  ■U-m +n  is  odd  or  if  sum  of  two 
indices  is  less  than  the  third, 


where  e=  |(^+m+n)  , A(k)  = 1*3  • • • (2k-l)/kl , k >0,  A(0)  = 1 . For 

programming  purposes  it  is  better  to  extend  the  definition  of  A(k)  by 

setting  A(k)  = 0,  k < 0,  in  which  case  the  two-index  sum  property  becomes 

a consequence  of  the  top  line  of  (28) . It  is  clear  that  P(£,m,n)  is 

symmetric  in  -C,m,n  . It  is  clear  that  P('t,m,n)  =0  if  -f+m+n  is  odd 

because  the  integrand  is  then  an  odd  function.  That  P(-t,m,n)  = 0 if, 

say,  t > m+n  follows  from  the  fact  that  the  Legendre  polynomial  P^  is 

orthogonal  to  all  polynomials  of  lower  degree,  in  particular,  P P 

m n 

It  is  important  that  (-t,m,n)  has  the  same  symmetry  and  vanishing 
properties  as  P(.l, m,n).  The  symmetry  is  obvious  from  the  definition 
(14a) . Again,  if  -t  > m+n  then 

- ter1 <p<-i  - pt,i>x  "5x dx 

by  (24) . But  the  factor  of  the  integrand  outside  the  parenthesis  is  a 
polynomial  of  degree  m+n-1,  and  both  P and  P^^  are  orthogonal  to 
it  since  -t+1  > -t^l  > m+n-1  . Finally,  the  same  formula  shows  that  if, 
say,  l is  odd  and  m and  n even  then  (-t,m,n)  = 0 since  the  inte- 
grand is  odd.  The  same  is  true  of  ■{.,«!, n all  odd.  This  exhausts  the 
cases  where  -f+m+n  is  odd  for  fixed  l in  tho  first  position.  By  sym- 
metry, then,  (-t,m,n)  = 0 for  -(+m+n  odd. 


Wo  now  wish  to  establish  tho  enumerative  formulas 
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(29) 


N3  + (N-2) 3 + 12(N2+(N-2)3)  + 20(N4N-2)1, 
i-T(N-l)3  + 12(N-1)2  + 20  (N-l)  1 


fox-  the  number  of  distinct  non-zero  (-t,m ,n)  with  -t,m,n  £ N with  N 
even  or  odd,  respectively,  i.e.,  in  view  of  the  symmetry,  the  number 
of  non-zero  (£,111, n)  with  N S £ ^ ra  ^ n.  We  reason  as  follows.  If 
N is  even,  then  there  are  N/2  (£,m,n),  L a m ^ n,  of  the  form  (N,N,2k), 

namely,  for  k = l,...,N/2;  N/2  of  the  form  (N,N-1 , 2k-l) , namely,  for 

k = l,...,N/2;  N/2  - 1 of  the  form  (N,N-2,2k),  namely,  for 
k = l,...,N/2-l;  N/2-2  of  the  form  (N,N-3,2k-l) , namely,  for 

k = 2, . . . , N/2-1}  etc.,  down  to  1 = N/2  - (N/2-1)  of  the  form  (N,N/2,N/2) 

for  a total  of  N/2  + (1  +...+  N/2)  = N/2  + (N/2) (N/2+l)/2  = N2/8+  3N/4. 


We  can  repeat  this  argument  with  -L  = N-2,  N-4,...,2  to  obtain 


i N 

8 Z n 
n even 


3 N 
| £ n 

n even 


. N/2 

- r 


„ N/2 

- £ 
2 . 
m=l 


= ^-(2  (N/2)  3 + 3 (N/2) 2 + N/2)  + N(N+2) 

= ^g(N34.  12N2  + 20N)  . 

Similarly,  there  are  (N-2)/2  of  the  form  (N-l, N-l, 2k) ; (N-2)/2  of 
the  form  (N-l, N-2, 2k-l) ; (N-2)/2  - 1 of  the  forni  (N-l,N-3,2k) , namely, 
for  k = l,...,(N-2)/2  - 1;  (N-2)/2  - 2 of  the  form  (N-l,N-4, 2k-l) , 
namely  for  k = 2, . . . , (N/2)/2  - 1;  and  so  on  down  to  1 = (N-2)/2  - 
(N/2  - 2)  of  the  form  (N-l, N/2, N/2-1)  for  a total  of  (N-2)/2  + 

(1  + ...  + (N-2)/2) . Repeating  the  argument  for  those  of  the  form 
(N-3,N-3,2k) , etc.,  one  sees  that  the  grand  total  for  -t  = N-l,N-3,  . . . ,3 
is  tho  same  as  that  obtained  above  but  with  N replaced  by  N-2.  Hence, 


for  N oven  tho  total  is  the  first  lino  of  (29).  But  a reexamination  of 
tho  last  part  of  tho  combinatorial  reasoning  shows  that  we  have  established 


I 
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tiiat  the  number  of  non-zoro  (-t,m ,n),  -l  2 m S n,  with  odd  *6  £ N,  N odd, 
is  precisely  equal  to  the  number  of  non-zero  (-t,m,n),  t a m 2 n,  with 
even  -t  *'  N-l.  Using  this  fact  we  see  that  for  N odd  we  get  a total 


of  twice  the  last  mentioned  number,  onco  for  N and  once  for  N-l,  or 
the  socond  line  of  (29) . 

To  completo  our  woi'k  it  is  only  necessary  to  derive  the  formula 
(19)  of  Section  3.  The  limits  in  all  the  integrals  in  (5b)  are  § = 0, 


= tt.  The  integral  in  the  first  term  is  thus  -2,  and  the  use  of  (17) 


accounts  for  the  first  term.  A typical  term  in  the  integral  of  3 cos  f 
in  the  second  term  of  (5b)  becomes,  when  one  sets  x = cos  § , 


, r1  n * 

-a.  J *4-  -s  — 


,n+l. 


n 


-1 


Jl-x* 


, . f n+1  n 1.,  , . . l.. 

dx  = -A  - — • + - — - (1  - (-1)  ) 

nL2n+l  2n+lJ 


= -(1  - (-l)n+1)A 


Here  we  have  used  the  consequence 


xdP  /dx  = [ (n+l)/(2n+l)ldP  /dx  + [ n/(2n+l)1 dP  /dx 
n n n+l 


of  (20)  and  (23) . For  the  next  term  of  (19)  we  take  the  next  two  terms 
of  (5b)  and  integrate  the  first  by  parts,  noting  that  the  integrated  out 


terms  disappear.  Then  the  same  device  as  used  just  now  and  (17)  complete 

2, 


the  term.  The  integral  of  3 sin  § cos  § in  the  last  linear  term 
becomes 

- £ A f P1(l  - x2)  (xA/i-x2)dx 

1 n 

1 


i n , n 
n=l  -1 


E A f P (1  - 3x  )dx  = (4/5)A 
. n a*  , n / 

n=l  -1 


since  1 - 3x  = 2P  , where  (7a)  and  an  integration  by  parts  have  been 
2 


used. 


I 


i 
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Tho  computation  in  (19)  of  the  first  two  non-linoar  terms  is  clear 

since  the  substitution  x = cos  £ turns  the  integrands,  with  the  aid 

* 

of  (17),  into  innor  products  of  3 with  itself  and  with  , respectively. 
The  same  dovices  convert  the  next  non-linear  integral,  except  for  the 
factor  K in  (18),  into 
(29) 


f — X-  dx  = 2 B*  A f x — x —?■  dx 

J -.2  _ _ ra  n J 


m,n=l 


1 dP 
x 

-1 


dP^ 
dx  x dx 


s * 1 (m-l+p(m-l))/2 


E B A f (raP  + 2 

..mnJ  m 41 
m,n=l  -1  i=l 


(41-l-2p(„-l))P21_2p(m.1)_1) 


(n-l+p(n-l))/2 

X (nP  + 2 

n . _ 

J=1 


(4j-2p(o-l)-l)P2J_2p(n.1)_1)  dx  , 


where  (26)  has  been  used,  which  gives  the  result  indicated  without  the 

addend  v in  parenthesis.  The  next  non-linear  term  in  (5b)  becomes 

(up  to  a constant  factor) 

(30) 

0 


r * / * 2 1 \ 

I t cos  § d§  = f ( £ B (-cot  ( P + m(m+l)cot  £ P j 
„ \ , m m m/ 

n tt  m=l 

X ( 2 A P1 ) sin  ? d§ 

\ n n / 
n=l 

1 / s . dP  dP  s 

= (’  ( 2 B A x2  ——  -p- ) dx 

J , \ , m n dx  dx/ 

-1  m,n=l 

1 / s dP  \ 

- f i 2 BA  m(m+l)  xP  ) dx  . 

J \ , m n m dx  / 

-1  m,n=l 

Here  the  first  equality  results  from  (9b),  (9c),  and  the  resulting  formula 

~ P1(cos  Q = -cot  § P1(cos  0 + m(m+l)P  (cos  9 • 
d § m m m 

But  tho  first  term  of  the  second  equality  in  (30)  is 
f ^*3  cot  § cos  § d§  , 


■ ■ 

• : 


7 . Discussion 

In  this  section  we  wish  to  take  up  a few  points  about  comparison 
with  different  theories,  other  results  and  possible  improvements  or  ex- 
tensions . 

First  of  nil,  in  using  E.  Reissner's  equations  in  small  finite 
deflection  form  and  retaining  only  the  terras  shown,  we  are  adopting  a 
quadrat ically  non-linear  theory  (quadratic  thoory,  for  shortness),  some- 
times called  a weakly  non-linear  theory.  While  this  is  admirably  adapted 
to  spectral  (Galerkin)  methods,  one  might  very  well  suspect  its  accuracy 
or  oven  relevance,  in  comparison  with  physical  reality  and  with  tho  full 


non-linear  equations  of  E.  Uoissner,  for  tho  largo  deflections  reported 
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3fC  * 

in  Section  5 in  tho  post  buckling  range  near  p = p and  on  the  rising 

Li 

section  of  tho  prossuro-deflcction  graph  to  tho  right.  While  conceding 
the  physical  unreality  of  the  rising  section  of  the  graph,  we  feel  thore 
are  some  grounds  for  a cautious  optimism  about  both  comparisons.  As  far 
as  the  fully  non-linear  Reissner  equations  go,  Mescall  [8]  has  reported 
only  insignificant  differences  between  the  predictions  of  that  theory  and 
the  quadratic  theory  for  large-deflection  buckling  of  spherical  caps,  and 
that  situation  and  ours  are  not  so  radically  different.  As  far  as  physical 
reality  goes,  from  the  very  beginnings  of  the  type  of  approximate  non-linear, 
in  fact,  quadratic  plate  and  shell  theories  we  are  using,  namely,  the 
Foppl-von  Karman  plate  equations,  comparisons  with  experiment  have  been 
remarkably  accurate  for  deflections  many  times  those  for  which  the  equa- 
tions were  originally  estimated  to  be  valid,  e.g.,  tens  of  thicknesses 
versus  one  thickness  for  the  F.-v.K.  equations.  Indeed,  there  is  almost 
a metamathematical  principle  of  unexpectedly  large  range  of  validity  in 
many  mathematical  models  in  physics  and  technology.  In  the  spherical  case 
the  shape  of  Figure  4 has  not  been  observed  experimentally  to  our  knowledge, 
but  that  may  well  be  duo  to  the  failure  of  the  elastic  stress-strain 
relations  over  portions  of  the  shell  at  snap-through.  In  view  of  the 
differential-geometric  plausibility  of  Figure  4 as  discussed  in  Section  5, 
one  can  conceive  of  a carefully  conducted  experiment,  possibly  with 
certain  constraints  during  buckling,  using  material  (not  metal)  of  ex- 
ceptional elastic  properties  (linear  elasticity  through  rather  large 


strains),  that  would  yield  the  shape  of  Figure  4. 

It  is  clear  that  tho  spectral  method  will  not  work  with  the  fully 
non-1 incar  Reissner  equations,  but  tho  pseudospectral  method  will.  The 
psoudospcctral  method  insorts  expansions  in  complete  function  sets  in  a 


set  of  differential  equations  but,  i*ather  than  reoxpanding  tho  equations, 


simply  demands  that  they  be  satisfied  at.  a prescribed  number  of  points. 

f i 

The  advantages  of  that,  from  file  computing  point  of  view  is  evident  oven 
for  the  cases  when  spectral  moLhods  are  applicable.  One  suspects, 
however,  a loss  of  accuracy  per  given  number  of  modes  vis-a-vis  spectral 
methods  in  those  cases. 

Our  computations  have  indicated  that  a/h  = 200  is  about  the  practical 
upper  limit  for  our  method  with  our  present  programs.  However,  it  would 
be  interesting  to  explore  the  possibility  of  applying  fast  transform 
methods  as  used  by  Orszag  [ 9]  in  conjunction  with  spectral  methods  in 
fluid  dynamics.  His  requirement  of  a finite  expansion  condition  is  precisely 
met  here.  The  resulting  shortening  of  the  computation  of  the  sums  in  (A') 
and  (B’>  might  well  make  possible  many  mode  computations,  necessary  for 
200  < a/h  £ 1000,  in  less  time  than  presently  required  for  a/h  < 100. 

Finally,  we  show  how  to  include  the  effect,  of  small  axisymmetric 
imperfections  in  Reissner’s  equations.  The  modified  equations  can  be 
solved  numerically  by  our  codes  after  making  the  necessary  obvious  modifi- 
cations. The  result  would  be  a modification  of  graphs  to  be  Figures  2 and 
3 near  the  origin  (clearly  one  would  need  expanded  deflection  scales  there).. 
Unfortunately  we  did  not  possess  tho  modified  equations  until  the  con- 
clusion of  the  numerical  work  reported  here,  and  we  do  not  have  numerical 
results  for  them.  Tho  important  thing  to  notice  is  that  the  imperfection 
sensitivity  (snap -through)  analysis  and  deep  post  buckling  analysis  appear 
simultaneously  as  features  of  one  non-linear  large  deflection  analysis. 

There  is  ample  precedent  for  this  observation  in  plate  and  shallow  shell 
analysis  in  the  literature,  hut  it  bears  reemphasizing. 

The  modified  equations  are  obtained  by  adding  tho  terms  PV  cot  £ “ 
a Dn^P  and  -P3  cot  £ to  the  right  sides  of  (4a)  and  (4b),  respectively, 


who  re 
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s 


£ e P (cos  f) 
n-1  n n 


and  e , 1 £ n - s,  is  the  (small)  imperfection  parameter  for  tho  n-th 

mode  imperfection.  These  terms  were  suggested  by  equations  (26-28)  of 

[ 2 1 . They  can  be  justified  by  inserting  cp^  = § - P in  equations  (III) 

and  (IV)  of  [ 12], getting  sin  = sin  ? - P cos  § , cos  = cos  § + 

2 

P sin  § there,  and  discarding  the  terms  P3  ^ , § P3  on  the  grounds  that 
they  are  essentially  cubic  and,  hence,  negligible  in  a quadratic  theory. 
The  result  is  to  add  the  terms 

'(£,  in , n)  ' 


K[ 


I z ritoll  e,  A , - 2 S 

n(n+l)  -1-v]  <£,,m=l  L Jn  J ^ m -t,m=l  L Jn 


e B*  + p%  , 
£ m n 


to  the  right  side  of  the  n-th  equation  of(A')and  loft  side  of  the  n-th 
equation  of  (B'),  respectively. 

The  following  simple  analysis  based  on  (A')  and  (B')  shows  the  im- 

* 

perfection-sensitivity  of  the  shell.  Take  s even,  A = B = 0 , 

n n 

n = 1, ... ,s-l  and  |a^ | so  small  that  its  cubes  are  negligible. 

* 

Eliminate  B to  obtain,  after  ignoring  cubes  and  dividing  by  A , 
s s 


(31) 


P - 1 = - 


Kf  s ( s+ 1 ) — 1 — vl  j 


This  shows  that  p will  decrease  linearly  from  1 for  Ag  small  and 

positivo,  i.e.,  inward  deflected  poles,  with  the  indicated 

If  s is  critical  for  a/h,  i.e.,  (12)  holds  approximately  then  (31) 
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ilc  line  lion  in  thicknesses,  to  not 


Thompson's  formula  f 13")  , 


P - 1 


•)  1 
r 3(i-\>  )i 
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